---
title: "Analysing data from study 1 telephone survey"
output: html_document
---
# 1. Basic system setup

```{r setup, include=FALSE}
knitr::opts_chunk$set(cache=FALSE)

# load libraries
library('tidyverse')
library('survey')
library('stargazer')
library('MASS')
library('knitr')
library('kableExtra')
library('psych')
library('nFactors')
library('RcmdrMisc')
library('moments')
library('ggeffects')
library('mediation')
library('ordinal')

# set ggplot theme
theme_set(theme_gray() +
            theme(panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            panel.background = element_rect(colour = "black",
                                            size = 1, fill = NA)))
```


# 2. Loading data: Telephone survey

```{r}
load(file = "replication_dataset_study1_telephone_poll.RData")
```

# 3. Check internal consistency for combined variables

## for grp_malleaby

Calculate the Cronbach's alpha with complete cases
```{r}
alpha_grp_malleaby <- telepoll %>% dplyr::select(V11:V14) %>% filter(complete.cases(.))
psych::alpha(alpha_grp_malleaby)$total$std.alpha
```

To calculate omega (an alternative measure of internal consistency)
```{r}
MBESS::ci.reliability(data=alpha_grp_malleaby, type="omega", conf.level = 0.95,
interval.type="bca", B=1000)
```


To check the number of factors for the components of grp_malleaby
### factor analysis

```{r}
n.factors <- 1

fit <- factanal(alpha_grp_malleaby, 
                n.factors,                # number of factors to extract
                scores=c("regression"),
                n.obs=nrow(alpha_grp_malleaby))

print(fit, digits=2, cutoff=.3, sort=TRUE)

```


```{r}
corMat <- cor(alpha_grp_malleaby)

(faPC  <- fa(r=corMat, nfactors=1, rotate="varimax", n.obs = nrow(alpha_grp_malleaby)))
```


```{r}

# Determine Number of Factors to Extract
ev <- eigen(cor(alpha_grp_malleaby)) # get eigenvalues
ap <- parallel(subject=nrow(alpha_grp_malleaby),
               var=ncol(alpha_grp_malleaby),
               rep=100,
               cent=.05)
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea)
plotnScree(nS)

```

## for ntl_homogty

Calculate the Cronbach's alpha with complete cases
```{r}
alp_ntl_homogty <- telepoll %>% dplyr::select(V10,V15,V18,V23) %>% filter(complete.cases(.))
psych::alpha(alp_ntl_homogty)$total$std.alpha
```

To calculate omega (an alternative measure of internal consistency)
```{r}
MBESS::ci.reliability(data=alp_ntl_homogty, type="omega", conf.level = 0.95,
interval.type="bca", B=1000)

```


## for stereotypes

Calculate the Cronbach's alpha with complete cases
```{r}
alp_stereotypes <- telepoll %>% dplyr::select(V20:V21) %>% filter(complete.cases(.))
psych::alpha(alp_stereotypes)$total$std.alpha

```

To calculate omega (an alternative measure of internal consistency)
```{r}
MBESS::ci.reliability(data=alp_stereotypes, type="omega", conf.level = 0.95,
interval.type="bca", B=1000)

```


# 4. Frequency distribution of attitudes towards migrant

## Appendix E: FREQUENCY DISTRIBUTION OF MIGRATION ATTITUDES IN STUDIES 1 

The distribution of attitudes towards migrants (V5) before recoding:
```{r}

ggplot(data = na.omit(telepoll), aes(x = as.numeric(V5))) +
  geom_histogram(binwidth = 0.5) +
  labs(x = "Attitudes towards migrants")

```

Skewness of the distribution of attitudes towards migrants before recoding:
```{r}
skewness(as.numeric(telepoll$V5), na.rm = TRUE)
```

Recode the attitudes towards migrants to binary
```{r}
# V5 is attitudes towards migrants
# combine 1-2; 3-5
telepoll$V5_bi <- telepoll$V5
telepoll$V5_bi[telepoll$V5_bi >= 1 & telepoll$V5_bi <= 2 & !is.na(telepoll$V5_bi)] <- 1
telepoll$V5_bi[telepoll$V5_bi >= 3 & telepoll$V5_bi <= 5 & !is.na(telepoll$V5_bi)] <- 2
telepoll$V5_bi <- factor(telepoll$V5_bi, levels = c(1:2))
```

# 5.. Summary statistics ------------------------------------------------------

## Appendix B: Summary statistics table.

a descriptive statistics for a list of variables in a table format.

```{r, results='asis'}
# create data frame for summary statistics
summary_stat <- as.data.frame(telepoll)
summary_stat <- summary_stat %>%
                        dplyr::select(V5, V5_bi, ntl_homogty, grp_malleaby,
                        stereotypes,
                        V16, V19, V6, V7, V22, V24,
                        V27, V29, V28, V26, V25,
                        w1, w2)
# convert factor into numeric
summary_stat$V5 <- summary_stat$V5 %>% as.numeric()
summary_stat$V5_bi <- summary_stat$V5_bi %>% as.numeric()
summary_stat$V28 <- summary_stat$V28 %>% as.numeric()
colnames(summary_stat) <- c("V5",
                            "attitu_immigrt",
                            "ntl_homogty",
                            "grp_malleaby",
                            "stereotypes",
                            "social_distance",
                            "cultural_threat",
                            "economic_threat",
                            "satisfied_hk_econ",
                            "proud_HKer",
                            "contact_immigrt",
                            "education",
                            "income",
                            "occupation",
                            "age",
                            "sex",
                            "w1",
                            "w2")

summary_stat_print <- summary_stat[, -c(1, length(summary_stat)-1, length(summary_stat))]

colnames(summary_stat_print) <- c("Attitudes towards migrants",
                            "National homogeneity",
                            "Group malleability",
                            "Stereotypes",
                            "Social distance",
                            "Cultural threat",
                            "Economic threat",
                            "Satisfaction with HK economy ",
                            "Proud of ��HKer�� identity",
                            "Contact with migrants",
                            "Education",
                            "Income",
                            "Occupation",
                            "Age",
                            "Sex")

stargazer(summary_stat_print, type = "text", summary = TRUE, rownames = FALSE)

```

# 6.. Correlations ------------------------------------------------------------


```{r}
corr_table <- t(cor(summary_stat_print,
                use = "complete.obs")[2:3, c(2,3,1,4:7)])
```

## Table 1: Bivariate Associations Between Perceived Group Malleability and National Homogeneity and Migration-Related Attitudes in Study 1

```{r, results='asis'}
stargazer(corr_table, type = "text", digits = 2)
```


## Table 1 (corresponding p-value for Table 1)

```{r}
cor_full <- rcorr.adjust(as.matrix(summary_stat_print[complete.cases(summary_stat_print),]),
                         type = c("pearson"),
                         use=c("complete.obs"))

cor_full$R$P[2:3, c(2,3,1,4:7)] %>% round(.,4) %>% t() %>% stargazer(type = "text", digits = 4)


```



# 7.. Regression analyses ---------------------------------

Create a data frame by selecting the variables needed, declare them either as ordered factor or factor.
```{r}
df.test <- summary_stat %>% 
  mutate_at(vars(-attitu_immigrt, -ntl_homogty, -grp_malleaby,
                 -stereotypes, -satisfied_hk_econ,
                 -proud_HKer, -age, -sex, -w1, -w2),
            as.ordered) %>% 
  mutate_at(vars(attitu_immigrt, sex), as.factor)
```

Run models

```{r}
m1.w <- lm(grp_malleaby ~
           ntl_homogty +
           satisfied_hk_econ +
           proud_HKer +
           contact_immigrt +
           education + 
           income +
           occupation +
           age +
           sex,
          weights = df.test$w2,
          data = df.test)

m2.w <- glm(attitu_immigrt ~
           ntl_homogty +
           satisfied_hk_econ +
           proud_HKer +
           contact_immigrt +
           education + 
           income +
           occupation +
           age +
           sex,
         weights = df.test$w2,
         data = df.test,
         family=binomial(link='probit'))

m3.w <- glm(attitu_immigrt ~
             grp_malleaby +
             satisfied_hk_econ +
             proud_HKer +
             contact_immigrt +
             education + 
             income +
             occupation +
             age +
             sex,
           weights = df.test$w2,
           data = df.test,
           family=binomial(link='probit'))

m4.w <- lm(stereotypes ~
             ntl_homogty +
             satisfied_hk_econ +
             proud_HKer +
             contact_immigrt +
             education + 
             income +
             occupation +
             age +
             sex,
           weights = df.test$w2,
           data = df.test)

m5.w <- lm(stereotypes ~
           grp_malleaby +
           satisfied_hk_econ +
           proud_HKer +
           contact_immigrt +
           education + 
           income +
           occupation +
           age +
           sex,
         weights = df.test$w2,
         data = df.test)

m6.w <- polr(social_distance ~
           ntl_homogty +
           satisfied_hk_econ +
           proud_HKer +
           contact_immigrt +
           education + 
           income +
           occupation +
           age +
           sex,
         weights = df.test$w2,
         data = df.test,
         method = "logistic",
         Hess = TRUE)

m7.w <- polr(social_distance ~
             grp_malleaby +
             satisfied_hk_econ +
             proud_HKer +
             contact_immigrt +
             education + 
             income +
             occupation +
             age +
             sex,
           weights = df.test$w2,
           data = df.test,
           method = "logistic",
           Hess = TRUE)

m8.w <- polr(cultural_threat ~
             ntl_homogty +
             satisfied_hk_econ +
             proud_HKer +
             contact_immigrt +
             education + 
             income +
             occupation +
             age +
             sex,
           weights = df.test$w2,
           data = df.test,
           method = "logistic",
           Hess = TRUE)

m9.w <- polr(cultural_threat ~
             grp_malleaby +
             satisfied_hk_econ +
             proud_HKer +
             contact_immigrt +
             education + 
             income +
             occupation +
             age +
             sex,
           weights = df.test$w2,
           data = df.test,
           method = "logistic",
           Hess = TRUE)

m10.w <- glm(economic_threat ~
             ntl_homogty +
             satisfied_hk_econ +
             proud_HKer +
             contact_immigrt +
             education + 
             income +
             occupation +
             age +
             sex,
           weights = df.test$w2,
           data = df.test,
           family=binomial(link='probit'))

m11.w <- glm(economic_threat ~
             grp_malleaby +
             satisfied_hk_econ +
             proud_HKer +
             contact_immigrt +
             education + 
             income +
             occupation +
             age +
             sex,
           weights = df.test$w2,
           data = df.test,
           family=binomial(link='probit'))


```

## Table 2: Effects of Perceived Group Malleability and National Homogeneity on Migration-Related Attitudes in Study 1 (weighted)

```{r results='asis', warning=FALSE}
stargazer(m1.w,m2.w,m3.w,m4.w,m5.w,
            type = "text",
          notes = "Entries are coefficients from various types of regressions; standard errors are in parentheses. All observations are weighted",
          notes.align = "l",
          dep.var.labels = c("Group malleability",
                            "Attitudes towards migrants",
                            "Stereotypes"),
          covariate.labels = c("Perceived homogeneity",
                               "Group malleability",
                               "Satisfied with HK economy",
                               "Proud of HKer identity",
                               "Social contact: 1-3",
                               "Social contact: 4 or above",
                               "Education: Secondary school",
                               "Education: A-Level",
                               "Education: Higher education (non-degree)",
                               "Education: Degree or above",
                               "Income: between HKD10,000 and HKD25,0000",
                               "Income: HKD25,000 or above",
                               "Occupation: Service and clerical support workers",
                               "Occupation: Laborers",
                               "Occupation: Students",
                               "Occupation: Housewives",
                               "Occupation: Others",
                               "Age",
                               "Female"),
          column.sep.width = "1pt")
```

## Table 3: Effects of Perceived Group Malleability and National Homogeneity on Migration-Related Attitudes in Study 1

```{r results='asis', warning=FALSE}
stargazer(m6.w,m7.w,m8.w,m9.w,m10.w,m11.w,
            type = "text",
          notes = "Entries are coefficients from various types of regressions; standard errors are in parentheses. All observations are weighted",
          notes.align = "l",
          dep.var.labels = c("Social distance",
                            "Cultural threats",
                            "Economic treats"),
          covariate.labels = c("Perceived homogeneity",
                               "Group malleability",
                               "Satisfied with HK economy",
                               "Proud of HKer identity",
                               "Social contact: 1-3",
                               "Social contact: 4 or above",
                               "Education: Secondary school",
                               "Education: A-Level",
                               "Education: Higher education (non-degree)",
                               "Education: Degree or above",
                               "Income: between HKD10,000 and HKD25,0000",
                               "Income: HKD25,000 or above",
                               "Occupation: Service and clerical support workers",
                               "Occupation: Laborers",
                               "Occupation: Students",
                               "Occupation: Housewives",
                               "Occupation: Others",
                               "Age",
                               "Female"),
          column.sep.width = "1pt")
```


# 8.. Predicted immigration attitudes --------------------------------------

## Figure 1: Predicted support of telephone survey respondents for migrants. Predicted values with 95% confidence intervals.

```{r}
m3.w_predict <- ggpredict(m3.w, terms = "grp_malleaby")

# generate chart
ggplot(m3.w_predict, aes(x = x)) +
  geom_line(aes(y = predicted)) +
  geom_line(aes(y=conf.low), linetype = "dashed") +
  geom_line(aes(y=conf.high), linetype = "dashed") +
  labs(y = "Predicted probabilities of supporting more migrants",
       x = "Perceived Group Malleability")

```


# 9. Indirect effect of implicit beliefs --------------------------------------

```{r mediation model setup}
df.mediate <- df.test %>%
             dplyr::select(attitu_immigrt,
             grp_malleaby,
             ntl_homogty,
             satisfied_hk_econ,
             proud_HKer,
             contact_immigrt,
             education, 
             income,
             occupation,
             age,
             sex,
             w1,
             w2) %>% na.omit()

model.m <- lm(grp_malleaby ~
             ntl_homogty +
             satisfied_hk_econ +
             proud_HKer +
             contact_immigrt +
             education + 
             income +
             occupation +
             age +
             sex,
           weights = df.mediate$w2,
           data = df.mediate)
summary(model.m)
```

```{r}
model.y <- glm(attitu_immigrt ~
             grp_malleaby +
             ntl_homogty +
             satisfied_hk_econ +
             proud_HKer +
             contact_immigrt +
             education + 
             income +
             occupation +
             age +
             sex,
           weights = df.mediate$w2,
           data = df.mediate,
           family=binomial(link='probit'))
summary(model.y)
```

```{r mediation model}
out.1 <- mediate(model.m, model.y, sims = 1000, boot = TRUE,
                 treat = "ntl_homogty", mediator = "grp_malleaby")

summary(out.1)
```

## Table 4: Mediation Effects of Implicit Beliefs About Group Malleability on Support Toward Migrants in Study 1

```{r}
cma <- summary(out.1)
acme <- c(cma$d.avg, cma$d.avg.ci)
direct_effect <- c(cma$z.avg, cma$z.avg.ci) 
total_effect <- c(cma$tau.coef, cma$tau.ci) 
round(rbind(acme, direct_effect, total_effect),4)

```


Visualisation of average causal mediated effect (ACME), average direct effect (ADE) and total effect.

```{r mediation effect plot}
plot(out.1)
```


```{r}
sens.out <- medsens(out.1, rho.by = 0.1, effect.type = "indirect", sims = 1000)
summary(sens.out)
```

## Appendix H: Sensitivity analysis for study 1 (fig. S1)

```{r}
plot(sens.out)
```

```{r}

```

